Modeling transport of extended interacting objects with drop-off phenomenon

We study a deterministic framework for important cellular transport phenomena involving a large number of interacting molecules called the excluded flow of extended interacting objects with drop-off effect (EFEIOD). This model incorporates many realistic features of biological transport process including the length of biological “particles” and the fact that they can detach along the biological ‘tracks’. The flow between the consecutive sites is unidirectional and is described by a “soft” simple exclusion principle and by repelling or attracting forces between neighboring particles. We show that the model admits a unique steady-state. Furthermore, if the parameters are periodic with common period T, then the steady-state profile converge to a unique periodic solution of period T. Simulations of the EFEIOD demonstrate several non-trivial effects of the interactions on the system steady-state profile. For example, detachment rates may help in increasing the steady-state flow by alleviating traffic jams that can exist due to several reasons like bottleneck rate or interactive forces between the particles. We also analyze the special case of our model, when there are no forces exerted by neighboring particles, and called it as the ribosome flow model of extended objects with drop-off effect (RFMEOD), and study the sensitivity of its steady-state to variations in the parameters.


Introduction
There are many important biological transport phenomena where the driving force for the movement of particles depends upon the constant source of energy [1]. One of the most known examples of such a system is intracellular transport carried out by motor proteins which are also known as biological molecular motors. They utilize the free energy that is produced during the chemical hydrolysis of adenosine triphosphate molecules (ATP) by converting it into mechanical work for their movement along cytoskeleton protein filaments [2]. These transport processes are modeled as the unidirectional or bidirectional flow of biological "particles" (RNA polymerases, ribosomes, motor proteins) along an ordered chain of sites (DNA, mRNA, microtubules).
Experimental investigations show that in many complex cellular processes such as mRNA translation by ribosomes or intracellular transport carried by motor proteins, the particles are larger than their step sizes and they usually function in large groups and interact with one another by binding and repelling actions based on the state of its neighboring particles [3]. Also, it has been seen that in many of these transport processes, the biological particles may get detached along the tracks. For example, kinesin-family motor proteins get detached from the microtubule after every power stroke or when their path is blocked [4,5]. Defects in kinesinlinked transport may disrupt the functioning of nerve cells and can cause many serious diseases [6]. The neuron-wide system requires intracellular transport of cargo throughout complex neuronal morphologies and its transport malfunction is one of the indications of some neuronal diseases like Alzheimer's [7]. Therefore, deriving mathematical models of these dynamical biological phenomena is important and crucial for understanding the collective behavior of the movement of particles and unraveling its biophysical aspects in the context of synthetic biology and biomedical applications. These transport phenomena are usually studied using several multi-particle lattice gas models [3]. The totally asymmetric simple exclusion process (TASEP) is an important model in non-equilibrium statistical physics that has been used to model molecular motors traffic and ribosome flow during mRNA translation [8][9][10]. In this stochastic model, particles hop unidirectionally with a step size of one along an ordered lattice of sites. The biological particles have volume and thus cannot overtake one other, i.e. as long as a site remain covered by a particle, it is inaccessible to the other particles thus obeying the simple exclusion principle. TASEP and its various versions have been introduced to model realistic observed features of interactions, extended objects and dissociations [11][12][13][14][15][16]. Unfortunately, rigorous analysis of TASEP is nontrivial and exact solutions exist in special simplified cases like the model with the homogeneous rates. Moreover detailed TASEP-type models are analyzed via various approximations and time-consuming Monte Carlo computer simulations. Therefore, a deterministic mathematical model called ribosome flow model (RFM) obtained via a mean-field approximation of TASEP has been introduced that is both amenable to mathematical analysis using tools from systems and control theory [17,18] and is easy to simulate. The analysis holds for any set of feasible parameter values. The RFM and its various extensions have been used to analyze many biological processes including positive feedbacks [19], the effect of competition for shared resources in translation [20], extended length of particles [29], bidirectional flow with Langmuir kinetics [21], networks of interconnected mRNAs [22], etc. Also, the analysis of the TASEP and its various extensions always provide approximate results which become accurate as number of sites n goes to infinity whereas analysis of the RFM and its various extensions hold true for every n.
An extension of the RFM called excluded flow with local repelling and binding model (EFRBM) was introduced to study the nearest-neighbor interactions between the motor proteins [23]. The advantage of this model is that it is that is amenable to rigorous analysis even in case of non-homogeneous rates. The EFRBM is a non-linear, continuous-time compartmental model for the unidirectional flow of interacting particles along a one-dimensional chain of n consecutive compartments. In this model, each particle is assumed to cover a single site and the nearest-neighbor effect is modeled by two interaction parameters q and r. The rate of movement of particle from site i to i + 1 depends upon the parameter q � 0 [r � 0] if site i + 2 [i − 1] is already occupied. A value q > 1 [q < 1] implies that the particle at site i will be strongly attracted [repelled] to the particle at site i + 2. This means that particle will move faster from site i to i + 1 creating a new bond with particle at site i + 2. Similarly a value r > 1 [r < 1] represents the detachment[attachment] force at site i by the neighboring particle at site i − 1. It has been shown that the trajectories of EFRBM evolve on the compact and convex set C n ≔ [0, 1] n .
In this paper, we extend the EFRBM to include the fact biological "particles" cover several sites and are susceptible to detach at various sites along the lattice. We refer to this model as the excluded flow of extended interacting objects with drop-off effect (EFEIOD). Using the theory of contractive dynamical systems, we prove that EFEIOD always converges to a steadystate. This steady-state depends on the length of the lattice n, the particle size ℓ, transition rates λ i s, detachment rates α i s and the interaction parameters q and r but not on the initial conditions. We also prove that it entrains to periodic excitations in the transition/detachment rates and the interaction parameters. This is important for the proper functioning of the biological processes that are excited by the periodic events. Analysis and simulations highlight the role of the effect of interactions on the steady-state flow. For example, in the case of strong attractions from the neighboring particle at site i − ℓ, the flow of particles from site i to i + 1 gets reduced, therefore an increase in detachment rate of particles at site i − ℓ leads to an easy steady-state flow. In the absence of interactions, we analyze the transport phenomena of mRNA translation with ribosome drop-off and called it RFMEOD. We also show using simulations that the RFMEOD correlates well with the TASEP with extended objects and including the drop-off phenomenon.
The EFEIOD presented here is more general than the EFRBM as it includes biologically observed more features such as particles with extended length and phenomena of dissociation of particles along the tracks.
The remainder of the paper is organized as follows. Section 2 describes the mathematical model. The next section presents our main theoretical results and the effects of the nearest-neighbor interactions on the steady-state behavior. Section 4 describes the application of the EFEIOD to model mRNA translation with ribosome drop-off and allows us to understand how a change in one of the parameters affects protein production. The final section concludes and summarizes the paper. To increase the readability of the paper, all the proofs are placed in the Appendix.

Model
The EFEIOD is a nonlinear, continuous time, compartmental model for the unidirectional flow of biological "particles" of size ℓ directed from left to right on a one-dimensional chain of n consecutive compartments or sites along the track.
3. r � 0, is the attachment/detachment force between any two existing consecutive particles.
4. q � 0, is the attachment/detachment force between any two new consecutive particles.
Each parameter λ i and α i has units of 1/time. A parameter q controls the repelling or binding forces between two new neighbors and a parameter r between two existing neighbors. In many studies, creating and breaking of bonds between the nearest neighbors has been viewed as opposite chemical reactions [24]. So, it is assumed that q , where E denotes the interaction energy, by applying the detailed balance arguments. The position of the particle along the lattice is denoted by the site covered by the leftmost end of it and this part is referred as the reader. Thus, 'the reader is at site i' means that the particle is located at site i and covers the sites i, i + 1, . . ., i + ℓ − 1. Let x i (t) 2 [0, 1] denote the normalized reader density of the biological particle at site i at time t and let y i (t) 2 [0, 1] denote its normalized coverage density at site i at time t, i.e., The term 'normalized' here means that each x i (t) and each y i (t) takes value in the interval [0, 1] for all t � 0. The value zero [one] corresponds to completely empty [full]. The schematic explanation of a particle with size ℓ on the lattice is shown in Fig 1. Eq (1) implies that the total particle coverage at any site i is the summation of the reader densities of ℓ consecutive sites left to site i. The state variables x i (t) and y i (t) can be interpreted as the probability that site i is occupied and covered respectively at time t. Hence, x i and y i are dimensionless.   are no readers at sites i − ℓ, i + ℓ and i + ℓ + 1, the transition rate is λ i and detachment rate is α i . Upper-right: When there is a reader at site i + ℓ + 1 and site i − ℓ does not have, the transition rate is λ i q and detachment rate is α i . Middle-left: When there is reader at site i − ℓ and no readers at sites i + ℓ and i + ℓ + 1, the transition rate is λ i r and detachment rate is α i r. Middle-right: When there are readers at sites i − ℓ and i + ℓ, detachment rate is α i r 2 . Lower part: When there are readers at sites i − ℓ and i + ℓ + 1, the transition rate is λ i qr and detachment rate is α i r.
To state the dynamical equations describing the EFEIOD, we introduce more notation for simplicity. Let The dynamics of the EFEIOD is described by n nonlinear first order ordinary differential equations: where and Eq (4) implies that the change in the reader density at site i is the inflow f i−1 (x) from site i − 1 to site i minus the outflow f i (x) to site i + 1 minus the outflow g i (x) to the cell environment.
Eq (6) can be explained as follows. The term x i represents that the reader flow from site i to site i + 1 increases with the reader density at site i. The term (1 − w i+ℓ ) represents a "soft" version of the simple exclusion principle which implies that the flow increases with the 'vacancy' level at site i + ℓ i.e. as the density in any of the ℓ consecutive sites increases the reader flow from site i to site i + 1 gradually decreases. The term (1 + (q − 1)z i+ℓ+1 ) represents that the reader flow from site i to site i + 1 also depends upon the reader density at site i + ℓ + 1 and increases [decreases] if q > 1 [q < 1]. The particle at site i + ℓ + 1 will attract [q > 1] or repel [q < 1] the particle that move from site i to i + 1. Similarly, the term (1+ (r − 1)z i−ℓ ) represents that the flow into site i + 1 also depends upon the reader density at site i − ℓ.
The term g i (x) in Eq (7) represents the detachment of particles from the site i to the cell environment. If r > 1 [r < 1] then the particles at sites i − ℓ and i + ℓ repel [attract] the particle at site i and increases [decreases] its detachment from site i.
The output rate from site n at any time t is given by Note that in the particular case, r = 1, q = 1, ℓ = 1 and α i = 0, the model gets reduced to the RFM [17]. Clearly, in the case when the length of the biological particle is equal to the lattice length, i.e. ℓ = n, there are no role of interaction forces in the system. The splitting of interaction energy E between the creation and breaking processes is not unique. Like in [15], we also assume that E is equally split between the rates r and q.
Note that Eq (9) implies r = 1/q and have a simple physical meaning. If E > 0 then there are attractive interactions in the system, i.e. the particle moves faster creating a new pair [q > 1] and the process of breaking out of the pair is slowed down [r < 1]. Similarly, E < 0 implies that there are repulsive interactions in the system. The case E = 0 corresponds to the fact that there are no interactions in the system and then we have q = r = 1.
The next section analyzes the EFEIOD using tools from systems and control theory and in particular contraction theory.
The set C n is not an invariant set of the EFEIOD as shown in the above example. Therefore, it is relevant to define a state space which is an invariant set of the dynamics. We assume that any initial condition belongs to the state-space defined as: where P is the lower triangular square matrix of size n with all entries zero, except for the entries on the main diagonal and (ℓ − 1) diagonals below the main diagonal that are ones. Note that the set H is a compact and convex set. We have shown in the following subsections that H is a relevant state-space for EFEIOD. Let int(H) and @H denote the interior and boundary of H respectively. Let x(t, a) denote the solution of Eq (4) at time t � 0 for the initial condition a 2 H.

Invariance and persistence
The following result shows that H is an invariant set for the dynamics of the EFEIOD.

Proposition 1 If a 2 H then the solution of EFEIOD satisfies x(t, a) 2 H for all t � 0. For any a 2 @H, x(t, a) 2 int(H) for t > 0.
This implies that the trajectories that emanate from the boundary of H 'immediately' enter the interior of H. The next proposition is useful because it shows that the solutions of the EFEIOD get 'immediately' uniformly separated from the boundary of H.

Proposition 2. For any τ > 0, there exists a compact and convex set H τ that is strictly contained in H such that for any a 2 H, x(t, a) 2 H τ , for all t � τ.
This means in particular that for any τ > 0 there exists d ¼ dðtÞ 2 0; 1

Contraction
Differential analysis provides a very useful way to study the behavior of certain non-linear dynamical systems. In particular, contraction theory is based on analyzing the time evolution of the distance between the trajectories that emanate from different initial conditions and have its applications to synchronization and reaction-diffusion partial differential equations [24,25]. However, for our proposed model, we needed a generalized version of contraction theory that has been defined in [26]. Consider a time-varying dynamical system, with the state x evolving on a compact and convex set O � R n . Let x(t, t 0 , a) denote the solution of (13) at time t � t 0 with x(t 0 ) = a. We say that (13) is contractive on O with respect a norm j:j : Now, the following generalization of contraction is required to apply contraction theory to the EFEIOD. The system (13) is said to be contractive after a small overshoot and short transient (SOST) [26] with respect a norm j:j : The next result shows that the EFEIOD satisfies this generalization of contraction. Let j:j 1 : R n ! R þ denote the L 1 norm, i.e. for x 2 R n , |x| 1 = |x 1 |+ |x 2 |+ . . .+ |x n |. Proposition 3. The EFEIOD is SOST with respect to the L 1 norm, i.e., for each � > 0 and each τ > 0 there exists c = c(τ, �) > 0 such that jxðt þ t; aÞ À xðt þ t; bÞj 1 � ð1 þ �Þ exp ðÀ tcÞja À bj 1 ; for all t � 0 and all a; b 2 H: ð16Þ This means that the EFEIOD is contractive after an arbitrarily small time transient τ and with an arbitrarily small overshoot (1 + �). This implies that any two initial feasible densities in the EFEIOD evolving in time become 'more similar' to each other at an exponential rate.

Global asymptotic stability
Since the convex and compact set H is an invariant set of the dynamics, it contains atleast one steady-state e [27]. By Proposition 1, we have e 2 int(H). Using Eq (16)  This means that, regardless of the initial density, all trajectories emanating from different initial conditions converge to the unique steady-state density that depends on the system parameters: transition rates λ i s, detachment rates α i s, interactions determined by q and r, particle size ℓ and length of the chain n. The next example demonstrates that the assumption q, r > 0 is necessary.
For q = r = 0 we have: Also, for q = 1 and r = 0, we have: Eqs (17) and (18) admits a continuum of steady-states, here [1 1 v] 0 is a steady-state for all v. Therefore, the assumption that q, r > 0 cannot be dropped.
The next example demonstrates the global asymptotic property, i.e. trajectories starting from different initial conditions in H asymptotically converge to a unique density profile along the lattice.
Example 3 Consider the EFEIOD with dimension n = 3, particle size ℓ = 2, rates λ i = 1, α i = 0.01, q = 1 and r = 1.  The next subsection analyze how the various parameters in the proposed model affects the steady-state output rate.

Analysis of the steady-state
At steady state, for x = e the left-hand side of all the equations in (4) is zero, so It follows from Eq (19) that if we multiply parameters λ i s and α i s by a scalar constant c > 0 then e will not change, i.e. e(cp) = e(p) where p = [λ 0 , λ 1 , . . ., λ n , α 1 , α 2 , . . ., α n ]. Also, R(cp) = cR(p), i.e. the output rate is homogeneous of order one w.r.t. the parameters λ i s and α i s. By Eq (19), we have: However, solving Eq (20), in general, is non-trivial. The next result shows that the derivatives of the equilibrium point coordinates with respect to the rates exists and are well defined. Let the mapping from the parameters to the unique equilibrium point be denoted by η, i.e., η i (γ) = e i , all i = 1, 2 . . ., n and γ = [λ 0 , λ 1 , . . ., λ n , α 1 , The above result allows to calculate the derivatives of the steady-state density if some of the parameters in the system are changed. This is useful to study the sensitivity of the steady-state w.r.t. small changes in the rates.

Effect of interactions
We demonstrate with several simulations the non-trivial effect of interactions on the steadystate of the EFEIOD.
The example below demonstrates that in the presence of strong attractive interactions, detachment of particles could be useful for increasing the flow of particles along the lattice.
Example 4 Consider the EFEIOD with dimension n = 9, particle size ℓ = 3, rates λ 0 = 1, λ i = 1 and α i = α. Fig 4 depicts that increasing the detachment rate increases the steady-state output for the higher values of q.
The above example suggests that for larger values of attractive interactions there could be a regulatory mechanism to increase the flow of particles in the system by allowing the particles to detach from the sites. The next example shows the positive role of increasing the detachment rate in the presence of a bottleneck rate at a site.
This can be explained as follows: for r < 1, a particle at site 5 will tend not to hop forward as there is strong attraction from particle at site 3. Therefore, allowing particles to detach from site 3, leads to an easy flow of particles from site 5 and this increases the flow. This is important to study as the interactions from the neighboring particles at the bottleneck rate further deteriorate the flow of particles along the lattice. In the case of no interactions i.e. q = 1 [r = 1], it can be seen that increasing the detachment rate leads to decrease in the steady-state flow which is always true as we theoretically analyze this special case in the next section.
The example above demonstrates that in the case of interactions, locally controlled detachment can avoid bottlenecks and can lead to faster movement of particles, and hence increasing the flow and alleviating the "traffic jams" [28]. One may perhaps think that increasing the particle size leads to a decrease in the steady-state output rate, but steady-state densities follow complicated behavior in the presence of interactions. It have been shown that when q = r = 1 and α i s = 0, steady-state output rate for ℓ > 1 is always less than steady-state output rate for the RFM [29]. But in the presence of interactions, increasing length does not always decrease the output rate as shown in the example below.
Example 6 Consider the EFEIOD with dimension n = 9, rates λ 0 = 1, λ i = 1 for all i, α i = 0. We vary the particle size ℓ. It can be seen in Fig 6 that for q = 13; for ℓ = 1, R = 0.0808, whereas for ℓ = 2, R = 0.1442 and for ℓ = 3, R = 0.1249. Furthermore, in the thermodynamical limit, i.e. number of sites goes to 1, the homogeneous case of TASEPEO with strong repulsions and particle size ℓ, and with entry and exit rates equal to one is in the maximal current phase, where the steady-state mean reader density 1=ð' þ 1 þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ' þ 1 p Þ and the steady-state output rate is 1=ðð1 þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ' þ 1 p Þ 2 Þ [14]. This implies that as ℓ goes to 1, the steady-state output and mean reader density go to zero. The next example shows that this is consistent with the results of our model. We define the steady-state mean reader density by r ¼ ð1=nÞ P n i¼1 e i . Example 7 Consider EFEIOD with dimension n = 100, rates λ 0 = 1, λ i = 1 for all i, α i = 0, q = 0.01, r = 100 and particle size ℓ. Fig 7 depicts that output rate R decrease with ℓ. Also, the steady-state mean reader density ρ decrease with ℓ as seen in Fig 8. The next example shows that in the presence of interactions, increase in an initiation rate does not always lead to increase in the output rate. However, in the case of no interactions i.e. q = 1 [r = 1], increase in an initiation rate due to feedbacks or due to an increase in number of 'free' biological particles leads to an increase in the steady-state output rate as we theoretically analyze this special case in the next section.
Example 8 Consider the EFEIOD with dimension n = 6, ℓ = 2, rates λ i = 1 for all i except λ 4 = 0.1 and α i = 0. We vary the initiation rate λ 0 . It can be seen in Fig 9 that steady-state output rate decreases with increase in λ 0 . Fig 10 depicts that the steady-state output rate increases with increase in λ 0 . Now, we analyze the effect of increasing the length of particle in the case q ! 1. Example 9 Consider the EFEIOD with dimension n = 3, λ 0 = 1, λ i = 1, α i = α, for all i. Fig  11 depicts that when q ! 1, the steady-state output rate decrease to zero. Fig 12 depicts that when q ! 1, the steady-state output rate saturates to a non-zero constant value depending on the value of α, i.e. R = 0.1910 for α = 0, R = 0.1720 for α = 0.5 and R = 0.1267 for α = 1.
A high value of q corresponds to a strong attachment between existing neighbors (small r) and a high tendency for creating new neighbors (large q), resulting in traffic jams and leading to a sharp decrease in the output rate. Therefore, the length of the particle has an interesting role to play in order to maintain a non-zero constant steady-state output rate in the case of weak repulsions.

Entrainment
There are many biological processes that are periodic [30], for example in translation-elongation mechanism; tRNA molecules [31], ATP levels [32], ribosome drop-off rate [33], translation initiation and elongation factors [34], oscillations in mRNA levels [35] and more may vary in a periodic manner and this results into the periodicity of the rates in the system. For the proper functioning of our body, certain biological systems must be in sync with the periodic changes induced due to the continuously changing environment [36,37]. Entrainment also plays an important part in designing extracellular biomedical systems [38]. An important question is: will the state variables of EFEIOD preserve the property of entrainment w.r.t. the parameters λ i s, α i s, q and r?
Assume that the λ i s, α i s, q and r are non-negative, uniformly bounded time-varying continuous functions satisfying: • There exists a (minimal) T > 0 such that every λ i s, α i s, q and r is a T-periodic function.

PLOS ONE
This model have been referred as periodic EFEIOD (PEFEIOD). The next result follows from the fact that EFEIOD is SOST on H and the known results on entrainment [39].

Theorem 2 The PEFEIOD admits a unique function � : R ! intðHÞ, that is T-periodic and for any initial condition a 2 H, the trajectory x(t, a) converges asymptotically to φ.
The above theorem implies that the state variables entrain to the periodic excitations in the parameters. The next example illustrates the behavior of PEFEIOD.
Example 10 Consider a PEFEIOD with dimension n = 3, ribosome size ℓ = 2, λ 0 = 1 λ i = 1, except for λ 2 (t) = 0.5 + 0.25 sin(πt/2) α i = 0.01, r = 5 and q = 1/5. Note, that there is single time-varying periodic rate in the network and all these rates are periodic with a common minimal period T = 4. We have taken two different initial conditions In general, to describe the effect of parameters on the system dynamics by a theoretical framework is cumbersome, as analyzing the set of non-linear equations that define the steadystate is not trivial. However, for a special case q = r = 1, the steady-state output rate sensitivity to variations in the parameters of the system can be answered rigorously. Moreover, the proposed general model was representative of the biology of molecular motors whereas this special case is important in the context of studying the ribosome flow and provides a tool in

Ribosome flow model with extended objects and ribosome dropoff
The synthesis of protein as directed by the mRNA template consisting of codons is carried out by ribosome and the process is referred to as translation [40]. The process broadly takes place in three steps: initiation where ribosomal complex assembles at the start codon of an mRNA chain; elongation where it moves along the mRNA in a forward series of steps forming a polypeptide chain of amino acids, and termination where it releases the chain that folds into functional protein and unbinds from the mRNA. Translation is a fundamental cellular process that occurs in all living beings at all times [41] and is known to consume most of the cell's energy [42]. Therefore, it is crucial to understand its dynamical aspects through mathematical modeling [43][44][45][46].
It is known from previous studies that the footprint of the ribosome on the mRNA is 10 to 20 codons [47][48][49]. Many ribosomes can simultaneously move on the same mRNA template, blocking the movement of other ribosomes behind it [50], resulting in traffic-like movement on the template, and these "traffic jams" are more severe in genes that are lowly expressed [51]. The ribosomes that initiate translation of mRNA sequence may not successfully complete it

PLOS ONE
and hence fails to produce a full-length protein product [52,53]. Hence, there are translational errors that can disrupt cellular fitness and can cause diseases [54]. Such errors can have multiple causes like ribosomal traffic jams, reading frameshifts [55], non-availability of tRNAs [56], misreading of codon, premature stop codons [57][58][59], etc. These errors often result in ribosome dissociating from the mRNA before reaching the stop codon called ribosome drop-off event, resulting in incomplete or incorrect peptides that are mostly non-functional, possibly toxic to the cell. The translational error due to premature translation termination seems to represent more than two-thirds of the overall errors and thus have a strong impact on protein formation [60,61]. Therefore, modeling mRNA translation with ribosome drop-off is important in analyzing the effect on the translation phenomena as it leads to a reduction in the rate of protein production.
To gain insights into these dynamical aspects of translation, we consider a special case of our model when q = r = 1 and we refer to this case as the ribosome flow model of extended objects with drop-off effect (RFMEOD). In this model, mRNA is treated as a one-dimensional lattice of length n, where n denotes the number of sites (codons) and every ribosome covers 1 � ℓ � n sites. The site 1 and n represents the start and stop codon respectively. The position of the ribosome along the mRNA is denoted by the site covered by the leftmost end of it. At any time t, if the leftmost edge of the ribosome is at site i, it means the reader is located at site i and the ribosome is translating site i and sites i, . . ., i + ℓ − 1 are covered by this ribosome. Ribosomes move unidirectionally from left to right by only one site on the template and no two ribosomes can occupy or cover the same site simultaneously.
The dynamics of RFMEOD is given by: . . .
. . . _ x n ¼ l nÀ 1 x nÀ 1 À l n x n À a n x n : The term λ i−1 x i−1 (1 − y i+ℓ−1 ) represents the reader flow from site i − 1 to site i. The flow increases with density level of readers at site i − 1 and decreases with coverage density y i+ℓ−1 = x i + x i+1 +. . .+ x i+ℓ−1 . The term α i x i represents the detachment of particles from the site i to the PLOS ONE cell environment. Also, the equations describing the last n − ℓ + 2 equations are linear, as a ribosome reading the last ℓ codons is the last particle and hence it move without any hindrance towards the last(stop) codon.
The output rate from site n at any time t, which is the protein production rate is given by RðtÞ ¼ ðl n þ a n Þx n ðtÞ: ð22Þ

Analysis of the steady-state
The RFMEOD Eq (21) can be written as where f 0 ðxÞ ≔l 0 ð1 À y ' Þ; f i ðxÞ ≔l i x i ð1 À y iþ' Þ; i ¼ 1; . . . ; n À 1; f n ðxÞ ≔l n x n ; Let R = (λ n + α n )e n denote the steady-state output rate. From Eq (25), we get This yields the following set of n + 1 equations in the n + 1 unknowns: e 1 , . . ., e n , R: and also Solving Eq (27) is in general non-trivial. Nevertheless, it can be solved in closed form in some special cases. Note that when α i = 0 for all i, RFMEOD gets reduced to RFMEO [29].

PLOS ONE
Example 11 Consider a RFMEOD with dimension n and with particle size ℓ = n. Consider homogeneous rates; λ i = λ, for i = 0, 1, . . ., n and α i = α, i = 1, 2, . . ., n. We have, and e i ¼ In case of totally homogeneous rates λ = α we have, Eq (26) can be used to prove various theoretical results. The next result shows that increasing any of the α i s, i 6 ¼ n, decreases R. In other words increasing any of the internal detachment rate decreases the steady-state protein production rate.
The next result shows that increasing any of the λ i s increases R. In particular, increasing the initiation rate always leads to an increase in protein synthesis rate. This result is consistent with a proposed canonical model of eukaryotic translation exhibiting a relation between initiation rate and protein expression [62].
Consider the case where all λ i = λ and α i = α. In this case, we can say more about steadystate densities.
Proposition 7 Consider a RFMEOD with dimension n and with particle size ℓ and λ i = λ and α i = α (6 ¼0). Then and This implies that the steady-state reader densities decrease between sites 1 and n and last ℓ sites reader density is given by Eq (31). The next example demonstrates this.
Example 12 The steady-state reader densities of the RFMEOD with dimension n = 16, for three particle sizes ℓ = 2, 4, 8 and λ 0 = 1, λ i = 1 and α i = 0.1, are depicted in Fig 14. It may be observed that steady-state reader densities monotonically decrease along the mRNA.
It have been seen that for fixed rates, the steady-state protein production rate in the RFMEO with ℓ > 1 is always less than the steady-state protein production rate with ℓ = 1 [29]. We also observed and investigated through simulations that this seems to hold true even in the case of the presence of drop-off phenomenon. The next example demonstrates that for fixed rates, steady-state output rate in the RFMEOD with ℓ > 1 is less than the steady-state output rate in the RFMEOD with ℓ = 1.
Example 13 Consider a RFMEOD with n = 300 sites, ribosome size ℓ and rates λ 0 = 0.8, λ i = 1, α i = 0.01, for all i. It can be seen that in Fig 15 that R monotonically decreases with ℓ.

RFMEOD with positive feedback
In eukaryotes, mRNA molecules sometimes form circular structures which promote recycling of the ribosomal subunits [63][64][65]. Therefore, it is biologically evident to include the fact that the translation initiation rate is affected by the premature and complete translation termination rate. This model can be used to fine tune the rate of protein production by ribosome-recycling in the case of changing ribosomal availability due to environmental stress [66]. We analyze the behavior of the RFMEOD as a control system after closing the loop from the output of ribosomes to input with positive linear feedback.
Here, the parameter k 1 represents the diffusion of ribosomes to the start codon of a mRNA molecule that is not related to recycling of ribosomes. The term k 2 ðl n x n þ P n i¼1 a i x i Þ represents the feedback due to recycling of ribosomes that have finished (partially or completely) the process of translating the mRNA as depicted in Fig 16. This is a generalization of the original RFMEOD as it includes both a term related to initiation rate with and without recycling of ribosomes. The next theorem proves that trajectories from any initial condition in H will always converge to a unique equilibrium point in H.

Theorem 3
The set H includes a unique equilibrium steady-state density e of the closed loop system (34). This equilibrium point is globally asymptotically stable in H, i.e. lim t!1 x(t, a) = e for any initial condition a 2 H. Example 14 Consider the closed loop system (34) with dimension n = 3, ℓ = 2, λ i = 1, α i = 0.01, k 1 = 1 and k 2 = 100. The next result provides information on the change of e w.r.t. the control parameters k 1 and k 2 .
Proposition 8 Suppose that the λ i s and α i s are fixed. Let e and � e correspond to the control parameters (k 1 , k 2 ) and ð � k 1 ; � k 2 Þ respectively. If k 1 ¼ � k 1 then e n < � e n if and only if k 2 < � k 2 . If k 2 ¼ � k 2 then e n < � e n if and only if k 1 < � k 1 . From a biophysical point of view, the above result inferred the intuitive result that increasing in any of the control parameters leads to an increase in the protein production rate. This result may be useful in context of biotechnology in order to improve levels of proteins in the host.

Validation through Monte Carlo simulations
Since the RFMEOD is a mean-field approximation of TASEP with extended objects and include an additional detachment rate at every site of the lattice, we ran MATLAB simulations of this process. A simulation begins with an empty chain of dimension n and continues for 10 8 time steps i.e. total simulation time. Each site can accommodate atmost one particle and a particle can only hop unidirectionally to a consecutive site if it is empty. The leftmost site that particle is covering is referred to as the reader. Every site i, i = 1, 2, . . ., n in the chain is associated with hopping rates λ i s and detachment rates α i s where the next hopping event time t k + � k or the next detachment event time t k + δ k is generated randomly. For site i, � k and δ k are random variables drawn from the exponential distribution with mean rate λ i and α i respectively. If hopping time is equal to the simulation time, then the reader at site i hops to site i + 1, provided site i + ℓ is empty. Similarly, if the detachment event time is equal to the simulation time then the reader dissociates from site i. The occupancy at each site is averaged throughout the simulations with the first 10 6 time steps discarded from the calculations to obtain the average steady-state reader density of each site.
In the example below, we show that simulations supports the modeling of dynamical aspects of translation using RFMEOD.

Discussion
In many biological processes like translation, cellular transport, gene transcription and many more, 'particles' move along one-dimensional "tracks". We studied a deterministic model called EFEIOD for the flow of particles along an ordered lattice of sites that encapsulates important cellular properties like detachment of particles from any site, nearest-neighbor interactions and the fact that most particles cover more than one site along the lattice. We analyzed this model using tools from systems and control theory, in particular contraction theory.
We proved that the EFEIOD converges to a unique steady-state density for any set of feasible parameters. In other words, EFEIOD is robust to the initial conditions. Moreover, we prove that if one or more of the parameters are time-varying periodic functions with a common period T, then the steady-state densities also converge to a periodic solution with period T. We demonstrate through simulations of the EFEIOD several useful observations. For example, increasing the particle size may sometimes lead to an increase in the output rate in the presence of weak repulsions. Surprisingly, we also show that increasing the detachment rate does not always decrease the output rate as elucidated in [28]. It is also important to note that several known models like RFM with positive feedback [19], RFMEO [29] and the model used in [67] for mRNA translation are special cases of the proposed model.
We also rigorously analyze a special case of the EFEIOD when q = r = 1 and called it RFMEOD to analyze the effects of ribosome drop-off on the translation process. The ribosome drop-off is important to study as it could significantly deteriorate the fitness of the host. We proved that increasing any one of the transition (detachment) rates of the RFMEOD always increases (decreases) the steady-state protein production and that in the homogeneous case, i.e. when all the transition rates are equal and all the detachment rates are equal, the reader density monotonically decreases along the lattice. We also modeled the observed phenomenon that many eukaryotic ribosomes may translate mRNA in multiples by including positive linear feedback in the RFMEOD.
The results reported can shed light on many biophysical properties of intracellular transport and may prove useful for applications in synthetic biology. One may consider to integrate another realistic feature of the cellular transport such as attachment of biological particles at different sites along the tracks in our model. Another research topic is studying the networks of EFEIOD and considering various phenomena like competition of resources in the network. We believe that the EFEIOD can be generalized to model and analyze more natural and artificial processes. Examples include coordination of large groups of organisms, traffic control and more.

Appendix
Proof of Proposition 1 and 2: The fact that H is an invariant set of the dynamics and have a repelling boundary follows from the equations from the EFEIOD. Let and with the z i s defined in Eq (2). Therefore, EFEIOD can be written as and with the w i s defined in Eq (3). Note that for r, q > 0 all the time-varying transition rates η i (t) are uniformly separated from zero and uniformly bounded and all the time-varying detachment rates are non-negative and uniformly bounded. Now, the proof of proposition follows from the results in [29,68].
Proof of Proposition 3: Let, and c n ðtÞ≔Z n ðtÞ ð1 À w nþ' ðtÞÞ: Now, combining the representations in Eqs 39 and 40 with the Eqs 37 and 38, we get Proposition 1 and the equations above imply that that EFEIOD can be interpreted as timevarying MFALK system with no backward and attachment dynamics with the well defined rates for all t > 0. Write the time-varying MFALK as _ x ¼ f ðx; tÞ with transition rates ψ i (t) and detachment rates δ i (t). A calculation shows that the Jacobian of f is J(t, x) = L(t, x) + D(t), where L is the matrix, and D is the diagonal matrix . . . ; À d nÀ 1 ; À c n À d n Þ: ð44Þ Hence Proposition 2 and the results in [21,68] imply that the EFEIOD is SOST on H and this completes the proof.
Proof of Proposition 4: It follows from the results in [21] and the argument used in the proof of proposition 3 in [23].
Now the proof follows by theorem 1.

Proof of Proposition 8:
We have equations for the RFMEOD with feedback at steady-state as follows: e n ≔ R l n ; e i ≔ R þ P nÀ 1 k¼iþ1 g k ðe k Þ l i ð1 À y iþ' Þ ; i ¼ n À 1; � � � ; 1; and also y ' ≔1 À R þ P nÀ 1 k¼1 g k ðe k Þ l 0 ðk 1 þ k 2 ðl n e n þ P n i¼1 a i e i ÞÞ : Suppose that k 1 ¼ � k 1 and k 2 < � k 2 . We have to prove that e n < � e n . We shall prove it by contradiction. Assume � e n � e n ð74Þ which implies that which further implies from Eq (73) that � e i � e i for all i ¼ 1; 2; � � � ; n À 1: Therefore, From Eq (73) and simplifying calculations we have, y ' À � y ' ¼ k 1 ðl n þ a n Þ ð � e n À e n Þ þ ðk 2 À � k 2 Þ ðl n ðl n þ a n Þ � e n þ a n ðl n þ a n Þ e n � e n Þ þðk 2 À � k 2 Þ ððl n þ a n Þ The fact that k 2 < � k 2 and Eqs (74) and (78) implies that which is a contradiction to Eq (77) and hence e n < � e n and the other part follows the same arguments.